#################################
### Revisiting White Backlash ###
###     Replication File      ###
###        AMT Data           ###
#################################

library(AER)
library(MASS)
library(multcomp)
library(stargazer)
library(texreg)

# Load clean data
load('AMTData.RData')

# Subset white respondents
WhiteData <- CleanData[CleanData$White == 1, ]

##############################
### Descriptive Statistics ###
###         pg. 7          ###
##############################

### --- Race --- ###
print(length(which(CleanData$Race == 1))) # white
#[1] 1653

print(length(which(CleanData$Race == 2))) # black
#[1] 103

print(length(which(CleanData$Race != 1 & CleanData$Race != 2))) # other
#[1] 375

### --- Gender --- ###
length(which(WhiteData$Gender == 0))/length(WhiteData$Gender) # female
#[1] 0.4679952

### --- Age --- ###
summary(as.factor(WhiteData$Age))
for(i in 1:7){
  print(length(which(as.factor(WhiteData$Age) == i))/length(as.factor(WhiteData$Age)))
}

# [1] 0.3683575
# [1] 0.321256
# [1] 0.1370773
# [1] 0.115942
# [1] 0.05012077
# [1] 0.004830918
# [1] 0.0006038647

### --- Education --- ###
# Education: 1=No HS, 2=HS, 3=Some college, 4=Associate degree, 5=Bachelors degree, 6=Masters degree, 7=Doctoral degree

for(i in 1:7){
  print(length(which(WhiteData$Educ == i))/length(WhiteData$Educ))
}

# [1] 0.003623188
# [1] 0.09359903
# [1] 0.2542271
# [1] 0.1038647
# [1] 0.3900966
# [1] 0.1189614
# [1] 0.03321256

### --- Income --- ###
#Income 1=<20K, 2=20-35, 3=35-50, 4=51-75, 5=75-100, 6=100+, 7=NA
for(i in 1:6){
  print(length(which(WhiteData$Income == i))/length(WhiteData$Income))
}

# [1] 0.1443237
# [1] 0.1944444
# [1] 0.1708937
# [1] 0.2071256
# [1] 0.1394928
# [1] 0.1256039

### --- Ideology --- ###
# Ideology: Conservative -> Liberal
for(i in 1:7){
  print(length(which(WhiteData$Ideol == i))/length(WhiteData$Ideol))
}

# [1] 0.06038647 - strong cons.
# [1] 0.1183575
# [1] 0.1123188
# [1] 0.1932367
# [1] 0.1449275
# [1] 0.2016908
# [1] 0.1672705 - strong lib.

### --- PID --- ###
# Rep -> Dem
# PID 7: Others coded in the Independent/Category
for(i in 1:7){
  print(length(which(WhiteData$PID7 == i))/length(WhiteData$PID7))
}

# [1] 0.08756039 - strong R
# [1] 0.1322464
# [1] 0.0718599
# [1] 0.1521739
# [1] 0.1437198
# [1] 0.2041063
# [1] 0.205314 - strong D

###################################
### Ordered Probit: Main Result ###
###        Table 3, p. 9        ###
###################################

# Treat1: 0 = none; 1 = Black; 2 = White
# Treat2: 0 = none; 1 = race frame
OP <- polr(as.factor(SupportDP) ~ Treat1*Treat2, data = WhiteData,
           method = 'probit', Hess = T)
# results
coeftest(OP)

# for table
texreg(OP, stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Black defendant',
                             'White defendant',
                             'Race frame',
                             'Black defendant $times$ race frame',
                             'White defendant $times$ race frame',
                             'Strongly Oppose $|$ Somewhat Oppose', 
                             'Somewhat Oppose $|$ Somewhat Favor', 
                             'Somewhat Favor $|$ Strongly Favor'))

# marginal effect of race frame in black photo condition
MEBlackPhoto <- glht(OP, linfct = c("Treat21 + Treat11:Treat21 = 0"))
summary(MEBlackPhoto)

# marginal effect of race frame in white photo condition
MEWhitePhoto <- glht(OP, linfct = c("Treat21 + Treat12:Treat21 = 0"))
summary(MEWhitePhoto)

##############################
### Bar Chart of % Support ###
###     Figure 1, p.10     ###
##############################

# add column of binary death penalty support
WhiteData$Favor <- ifelse(WhiteData$SupportDP <=2, 0, 1)

# generate means and confidence intervals for each treatment
NPNF <-  t.test(WhiteData[WhiteData$TreatmentFactor == 'NP_NF', 'Favor'])
BPNF <-  t.test(WhiteData[WhiteData$TreatmentFactor == 'BP_NF', 'Favor'])
WPNF <-  t.test(WhiteData[WhiteData$TreatmentFactor == 'WP_NF', 'Favor'])
NPBF <-  t.test(WhiteData[WhiteData$TreatmentFactor == 'NP_BF', 'Favor'])
BPBF <-  t.test(WhiteData[WhiteData$TreatmentFactor == 'BP_BF', 'Favor'])
WPBF <-  t.test(WhiteData[WhiteData$TreatmentFactor == 'WP_BF', 'Favor'])

# store vector of point estimates
PointEstimates <- c(NPNF$estimate, NPBF$estimate, BPNF$estimate,
                    BPBF$estimate, WPNF$estimate, WPBF$estimate)*100

# store vector of confidence intervals
ConfInts <- list(NPNF$conf.int, NPBF$conf.int, BPNF$conf.int, 
                 BPBF$conf.int, WPNF$conf.int, WPBF$conf.int)

# create barplot
barplot(PointEstimates, ylim = c(0, 70), space = c(0,0,1,0,1,0),
        xlab = '', ylab = 'Percent favoring death penalty', axes = F, 
        names.arg = F, col = c('white', 'grey65'), 
        main = '')
# add y-axis
axis(2, las = 1)
# add x-axis
axis(1, at = c(1, 4, 7), tick = F, line = -5, 
     labels = c('No Defendant', 'Black Defendant', 'White Defendant'), outer = T)
# add legend
legend('top', bty = 'n', title = 'Frame', cex = 0.65,
       legend = c('None', 'Race'), fill = c('white', 'grey65'))
# add confidence bars
for (i in 1:6) {
  segments(x0 = c(0, 1, 3, 4, 6, 7)[i]+0.5, x1 = c(0, 1, 3, 4, 6, 7)[i]+0.5, 
           y0 = ConfInts[[i]][1]*100, y1 = ConfInts[[i]][2]*100)
  segments(x0 = c(0, 1, 3, 4, 6, 7)[i]+0.25, x1 = c(0, 1, 3, 4, 6, 7)[i]+0.75,
           y0 = ConfInts[[i]][1]*100, y1 = ConfInts[[i]][1]*100)
  segments(x0 = c(0, 1, 3, 4, 6, 7)[i]+0.25, x1 = c(0, 1, 3, 4, 6, 7)[i]+0.75,
           y0 = ConfInts[[i]][2]*100, y1 = ConfInts[[i]][2]*100)
}
# add text
text(0.5,5, 'N=')
text(0.5,2, NPNF$parameter+1)
text(1.5,2, NPBF$parameter+1)
text(3.5,2, BPNF$parameter+1)
text(4.5,2, BPBF$parameter+1)
text(6.5,2, WPNF$parameter+1)
text(7.5,2, WPBF$parameter+1)

######################
### OLS: All Races ###
###   Table S1     ###
######################

# Treat1: 0 = none; 1 = Black; 2 = White
# Treat2: 0 = none; 1 = race frame

# White respondent OLS
OLS <- lm(SupportDP ~ Treat1*Treat2, data = WhiteData)

# results
summary(OLS)

# robust SEs
sqrt(diag(vcovHC(OLS)))

# subset black respondents 
BlackData <- CleanData[CleanData$Race == 2, ]

# Black respondent OLS
OLSBlack <- lm(SupportDP ~ Treat1*Treat2, data = BlackData)

# results
summary(OLSBlack)

# robust SEs
sqrt(diag(vcovHC(OLSBlack)))

# subset other respondents
OtherData <- CleanData[CleanData$Race != 1 & CleanData$Race != 2, ]

# other respondents OLS
OLSOther <- lm(SupportDP ~ Treat1*Treat2, data = OtherData)

# results
summary(OLSOther)

# robust SEs
sqrt(diag(vcovHC(OLSOther)))

# table output
texreg(list('Whites' = OLS, 
            'Blacks' = OLSBlack, 
            'Others' = OLSOther), 
       custom.coef.map = list('Treat21' = 'Race frame',
                              'Treat11' = 'Black defendant',
                              'Treat12' = 'White defendant',
                              'Treat11:Treat21' = 'Race frame $times$ black defendant',
                              'Treat12:Treat21' = 'Race frame $times$ white defendant',
                              '(Intercept)' = 'Constant'),
       override.se = list(sqrt(diag(vcovHC(OLS))), 
                       sqrt(diag(vcovHC(OLSBlack))),
                       sqrt(diag(vcovHC(OLSOther)))),
       stars = numeric(0), digits = 3)

####################################
### Bar Chart of Ordinal Support ###
###           Figure S1          ###
####################################

# generate means and confidence intervals for each treatment
NPNFOrd<-  t.test(WhiteData[WhiteData$TreatmentFactor == 'NP_NF', 'SupportDP'])
BPNFOrd <-  t.test(WhiteData[WhiteData$TreatmentFactor == 'BP_NF', 'SupportDP'])
WPNFOrd <-  t.test(WhiteData[WhiteData$TreatmentFactor == 'WP_NF', 'SupportDP'])
NPBFOrd <-  t.test(WhiteData[WhiteData$TreatmentFactor == 'NP_BF', 'SupportDP'])
BPBFOrd <-  t.test(WhiteData[WhiteData$TreatmentFactor == 'BP_BF', 'SupportDP'])
WPBFOrd <-  t.test(WhiteData[WhiteData$TreatmentFactor == 'WP_BF', 'SupportDP'])

# store vector of point estimates
PointEstimatesOrd <- c(NPNFOrd$estimate, NPBFOrd$estimate, BPNFOrd$estimate,
                       BPBFOrd$estimate, WPNFOrd$estimate, WPBFOrd$estimate)

# vector of confidence intervals
ConfIntsOrd <- list(NPNFOrd$conf.int, NPBFOrd$conf.int, BPNFOrd$conf.int, 
                    BPBFOrd$conf.int, WPNFOrd$conf.int, WPBFOrd$conf.int)

# create barplot
barplot(PointEstimatesOrd, ylim = c(1, 3), space = c(0,0,1,0,1,0), xpd = F,
        xlab = '', ylab = 'Average death penalty support', axes = F, names.arg = F, 
        col = c('white', 'grey65'), 
        main = '')
# add y-axis
axis(2, las = 1)
# add x-axis
axis(1, at = c(1, 4, 7), tick = F, line = -5, 
     labels = c('No Defendant', 'Black Defendant', 'White Defendant'), outer = T)
# add legend
legend('top', bty = 'n', title = 'Frame', cex = 0.65,
       legend = c('None', 'Race'), fill = c('white', 'grey65'))
# add confidence bars
for (i in 1:6) {
  segments(x0 = c(0, 1, 3, 4, 6, 7)[i]+0.5, x1 = c(0, 1, 3, 4, 6, 7)[i]+0.5, 
           y0 = ConfIntsOrd[[i]][1], y1 = ConfIntsOrd[[i]][2])
  segments(x0 = c(0, 1, 3, 4, 6, 7)[i]+0.25, x1 = c(0, 1, 3, 4, 6, 7)[i]+0.75,
           y0 = ConfIntsOrd[[i]][1], y1 = ConfIntsOrd[[i]][1])
  segments(x0 = c(0, 1, 3, 4, 6, 7)[i]+0.25, x1 = c(0, 1, 3, 4, 6, 7)[i]+0.75,
           y0 = ConfIntsOrd[[i]][2], y1 = ConfIntsOrd[[i]][2])
}
segments(x0=0, x1=2, y0=1, y1=1, lwd=2)
segments(x0=3, x1=5, y0=1, y1=1, lwd=2)
segments(x0=6, x1=8, y0=1, y1=1, lwd=2)
text(0.5,1.15, 'N=')
text(0.5,1.05, NPNFOrd$parameter+1)
text(1.5,1.05, NPBFOrd$parameter+1)
text(3.5,1.05, BPNFOrd$parameter+1)
text(4.5,1.05, BPBFOrd$parameter+1)
text(6.5,1.05, WPNFOrd$parameter+1)
text(7.5,1.05, WPBFOrd$parameter+1)

###################################
### Ordered Probit: PH Controls ###
###         Table S3            ###
###################################

# regression controlling for PH controls
OPControlled <- polr(as.factor(SupportDP) ~ BlackCrime*Treat1*Treat2
                     + GralCrime*Treat1*Treat2
                     + AntiBlack*Treat1*Treat2
                     + Fear*Treat1*Treat2
                     + Punit*Treat1*Treat2
                     + PID7*Treat1*Treat2
                     + Ideol*Treat1*Treat2
                     + Educ*Treat1*Treat2
                     + Gender*Treat1*Treat2
                     + Income*Treat1*Treat2
                     + Age*Treat1*Treat2,
                     data = WhiteData,
                     method = 'probit')

# results
coeftest(OPControlled)

texreg(OPControlled, stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Black defendant',
                             'White defendant',
                             'Race frame',
                             'Black crime attribution',
                             'General crime attribution',
                             'Anti-black stereotypes',
                             'Fear of crime',
                             'Punitiveness',
                             'Party ID',
                             'Ideology',
                             'Education',
                             'Female',
                             'Income',
                             'Age',
                             'Black defendant $\\times$ race frame',
                             'White defendant $\\times$ race frame',
                             'Strongly Oppose $|$ Somewhat Oppose', 
                             'Somewhat Oppose $|$ Somewhat Favor', 
                             'Somewhat Favor $|$ Strongly Favor'))

############################################
### Heterogeneous Treatment: Black Crime ###
###                Table S4              ###
############################################

# binary indicator for dispositional attribution of black crime
WhiteData$BCA <- WhiteData$BlackCrime > 2

# ordered probit conditioning on Black Crime Attribution
BCAOP <- polr(as.factor(SupportDP) ~ BCA*Treat1*Treat2,
              data = WhiteData, 
              method = 'probit')

# results                
coeftest(BCAOP)

# marginal effect of BCA on race frame only
MENoPhotoRaceFrame <- glht(BCAOP,
                           linfct = c("Treat21 + BCATRUE:Treat21 = 0"))
summary(MENoPhotoRaceFrame)
# power = 0.129
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$BCA == T & WhiteData$TreatmentFactor == 'NP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$BCA == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.01045/0.11797,
             sig.level = 0.05,
             alternative = 'two.sided')

# marginal effect of BCA on race frame for black defendant condition
MEBlackPhotoRaceFrame <- glht(BCAOP,
                              linfct = c("Treat21 + BCATRUE:Treat21 + Treat11:Treat21 + BCATRUE:Treat11:Treat21 = 0"))
summary(MEBlackPhotoRaceFrame)
# power = 1
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$BCA == T & WhiteData$TreatmentFactor == 'BP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$BCA == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.1441/0.1204,
             sig.level = 0.05,
             alternative = 'two.sided')

# marginal effect of BCA on race frame for white defendant condition
MEWhitePhotoRaceFrame <- glht(BCAOP,
                              linfct = c("Treat21 + BCATRUE:Treat21 + Treat12:Treat21 + BCATRUE:Treat12:Treat21 = 0"))
summary(MEWhitePhotoRaceFrame)
# power = 1
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$BCA == T & WhiteData$TreatmentFactor == 'WP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$BCA == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.1483/0.1169,
             sig.level = 0.05,
             alternative = 'two.sided')

# table output of ordered probit
texreg(BCAOP, stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Dispositional black crime attribution',
                             'Black defendant',
                             'White defendant',
                             'Race frame',
                             'Dispositional black crime attribution $times$ black defendant',
                             'Dispositional black crime attribution $times$ white defendant',
                             'Dispositional black crime attribution $times$ race frame',
                             'Black defendant $times$ race frame',
                             'White defendant $times$ race frame',
                             'Dispositional black crime attribution $times$ black defendant $times$ race frame',
                             'Dispositional black crime attribution $times$ white defendant $times$ race frame',
                             'Strongly Oppose $|$ Somewhat Oppose', 
                             'Somewhat Oppose $|$ Somewhat Favor', 
                             'Somewhat Favor $|$ Strongly Favor'))


############################################
### Heterogeneous Treatment: Stereotypes ###
###                Table S5              ###
############################################

# binary indicator for anti-black stereotypes
WhiteData$ABS <- WhiteData$AntiBlack > 0

# ordered probit conditioning on anti-black stereotypes
ABSOP <- polr(as.factor(SupportDP) ~ ABS*Treat1*Treat2,
              data = WhiteData, 
              method = 'probit')

# results                
coeftest(ABSOP)

# marginal effect of ABS on race frame only
MENoPhotoRaceFrame <- glht(ABSOP,
                           linfct = c("Treat21 + ABSTRUE:Treat21 = 0"))
summary(MENoPhotoRaceFrame)
# power = 0.850
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$ABS == T & WhiteData$TreatmentFactor == 'NP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$ABS == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.0901/0.1837,
             sig.level = 0.05,
             alternative = 'two.sided')

# marginal effect of ABS on race frame for black defendant condition
MEBlackPhotoRaceFrame <- glht(ABSOP,
                              linfct = c("Treat21 + ABSTRUE:Treat21 + Treat11:Treat21 + ABSTRUE:Treat11:Treat21 = 0"))
summary(MEBlackPhotoRaceFrame)
# power = 1
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$ABS == T & WhiteData$TreatmentFactor == 'BP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$ABS == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.2689/0.1775,
             sig.level = 0.05,
             alternative = 'two.sided')

# marginal effect of ABS on race frame for white defendant condition
MEWhitePhotoRaceFrame <- glht(ABSOP,
                              linfct = c("Treat21 + ABSTRUE:Treat21 + Treat12:Treat21 + ABSTRUE:Treat12:Treat21 = 0"))
summary(MEWhitePhotoRaceFrame)
# power = 0.999
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$ABS == T & WhiteData$TreatmentFactor == 'WP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$ABS == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.1825/0.1770,
             sig.level = 0.05,
             alternative = 'two.sided')

# table output of ordered probit
texreg(ABSOP, stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Anti-black stereotypes',
                             'Black defendant',
                             'White defendant',
                             'Race frame',
                             'Anti-black stereotypes $times$ black defendant',
                             'Anti-black stereotypes $times$ white defendant',
                             'Anti-black stereotypes $times$ race frame',
                             'Black defendant $times$ race frame',
                             'White defendant $times$ race frame',
                             'Anti-black stereotypes $times$ black defendant $times$ race frame',
                             'Anti-black stereotypes $times$ white defendant $times$ race frame',
                             'Strongly Oppose $|$ Somewhat Oppose', 
                             'Somewhat Oppose $|$ Somewhat Favor', 
                             'Somewhat Favor $|$ Strongly Favor'))


###########################################
### Heterogeneous Treatment: Republican ###
###                Table S6             ###
###########################################

# binary indicator for Republican party ID
WhiteData$GOP <- WhiteData$PID7 < 4

# ordered probit conditioning on Republican partisanship
GOPOP <- polr(as.factor(SupportDP) ~ GOP*Treat1*Treat2,
              data = WhiteData, 
              method = 'probit')

# results                
coeftest(GOPOP)

# marginal effect of GOP on race frame only
MENoPhotoRaceFrame <- glht(GOPOP,
                           linfct = c("Treat21 + GOPTRUE:Treat21 = 0"))
summary(MENoPhotoRaceFrame)
# power = 0.725
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$GOP == T & WhiteData$TreatmentFactor == 'NP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$GOP == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.06645/0.17007,
             sig.level = 0.05,
             alternative = 'two.sided')

# marginal effect of GOP on race frame for black defendant condition
MEBlackPhotoRaceFrame <- glht(GOPOP,
                              linfct = c("Treat21 + GOPTRUE:Treat21 + Treat11:Treat21 + GOPTRUE:Treat11:Treat21 = 0"))
summary(MEBlackPhotoRaceFrame)
# power = 0.085
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$GOP == T & WhiteData$TreatmentFactor == 'BP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$GOP == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.01686/0.18196,
             sig.level = 0.05,
             alternative = 'two.sided')

# marginal effect of GOP on race frame for white defendant condition
MEWhitePhotoRaceFrame <- glht(GOPOP,
                              linfct = c("Treat21 + GOPTRUE:Treat21 + Treat12:Treat21 + GOPTRUE:Treat12:Treat21 = 0"))
summary(MEWhitePhotoRaceFrame)
# power = 1
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$GOP == T & WhiteData$TreatmentFactor == 'WP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$GOP == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.2046/0.1689,
             sig.level = 0.05,
             alternative = 'two.sided')

# table output of ordered probit
texreg(GOPOP, stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Republican',
                             'Black defendant',
                             'White defendant',
                             'Race frame',
                             'Republican $times$ black defendant',
                             'Republican $times$ white defendant',
                             'Republican $times$ race frame',
                             'Black defendant $times$ race frame',
                             'White defendant $times$ race frame',
                             'Republican $times$ black defendant $times$ race frame',
                             'Republican $times$ white defendant $times$ race frame',
                             'Strongly Oppose $|$ Somewhat Oppose', 
                             'Somewhat Oppose $|$ Somewhat Favor', 
                             'Somewhat Favor $|$ Strongly Favor'))


#############################################
### Heterogeneous Treatment: Conservative ###
###                 Table S7              ###
#############################################

# binary indicator for ideological conservatism
WhiteData$Con <- WhiteData$Ideol < 4

# ordered probit conditioning on conservatism
ConOP <- polr(as.factor(SupportDP) ~ Con*Treat1*Treat2,
              data = WhiteData, 
              method = 'probit')

# results                
coeftest(ConOP)

# marginal effect of conservatism on race frame only
MENoPhotoRaceFrame <- glht(ConOP,
                           linfct = c("Treat21 + ConTRUE:Treat21 = 0"))
summary(MENoPhotoRaceFrame)
# power = 0.082
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Con == T & WhiteData$TreatmentFactor == 'NP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Con == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.0132/0.1677,
             sig.level = 0.05,
             alternative = 'two.sided')

# marginal effect of conservatism on race frame for black defendant condition
MEBlackPhotoRaceFrame <- glht(ConOP,
                              linfct = c("Treat21 + ConTRUE:Treat21 + Treat11:Treat21 + ConTRUE:Treat11:Treat21 = 0"))
summary(MEBlackPhotoRaceFrame)
# power = 0.951
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Con == T & WhiteData$TreatmentFactor == 'BP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Con == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.1055/0.1755,
             sig.level = 0.05,
             alternative = 'two.sided')

# marginal effect of conservatism on race frame for white defendant condition
MEWhitePhotoRaceFrame <- glht(ConOP,
                              linfct = c("Treat21 + ConTRUE:Treat21 + Treat12:Treat21 + ConTRUE:Treat12:Treat21 = 0"))
summary(MEWhitePhotoRaceFrame)
# power = 1
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Con == T & WhiteData$TreatmentFactor == 'WP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Con == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.2769/0.1786,
             sig.level = 0.05,
             alternative = 'two.sided')

# table output of ordered probit
texreg(ConOP, stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Conservative',
                             'Black defendant',
                             'White defendant',
                             'Race frame',
                             'Conservative $times$ black defendant',
                             'Conservative $times$ white defendant',
                             'Conservative $times$ race frame',
                             'Black defendant $times$ race frame',
                             'White defendant $times$ race frame',
                             'Conservative $times$ black defendant $times$ race frame',
                             'Conservative $times$ white defendant $times$ race frame',
                             'Strongly Oppose $|$ Somewhat Oppose', 
                             'Somewhat Oppose $|$ Somewhat Favor', 
                             'Somewhat Favor $|$ Strongly Favor'))


##############################################
### Heterogeneous Treatment: Low Education ###
###                 Table S8               ###
##############################################

# binary indicator for low education
WhiteData$LowEdu <- WhiteData$Educ < 4

# ordered probit conditioning on low education
LowEduOP <- polr(as.factor(SupportDP) ~ LowEdu*Treat1*Treat2,
                 data = WhiteData, 
                 method = 'probit')

# results                
coeftest(LowEduOP)

# marginal effect of low ed on race frame only
MENoPhotoRaceFrame <- glht(LowEduOP,
                           linfct = c("Treat21 + LowEduTRUE:Treat21 = 0"))
summary(MENoPhotoRaceFrame)
# power = 1
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$LowEdu == T & WhiteData$TreatmentFactor == 'NP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$LowEdu == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.1459/0.1459,
             sig.level = 0.05,
             alternative = 'two.sided')

# marginal effect of low ed on race frame for black defendant condition
MEBlackPhotoRaceFrame <- glht(LowEduOP,
                              linfct = c("Treat21 + LowEduTRUE:Treat21 + Treat11:Treat21 + LowEduTRUE:Treat11:Treat21 = 0"))
summary(MEBlackPhotoRaceFrame)
# power = 0.999
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Con == T & WhiteData$TreatmentFactor == 'BP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Con == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.1377/0.1567,
             sig.level = 0.05,
             alternative = 'two.sided')

# marginal effect of low ed on race frame for white defendant condition
MEWhitePhotoRaceFrame <- glht(LowEduOP,
                              linfct = c("Treat21 + LowEduTRUE:Treat21 + Treat12:Treat21 + LowEduTRUE:Treat12:Treat21 = 0"))
summary(MEWhitePhotoRaceFrame)
# power = 0.370
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Con == T & WhiteData$TreatmentFactor == 'WP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Con == T & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.04043/0.15957,
             sig.level = 0.05,
             alternative = 'two.sided')

# table output of ordered probit
texreg(LowEduOP, stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Low education',
                             'Black defendant',
                             'White defendant',
                             'Race frame',
                             'Low education $times$ black defendant',
                             'Low education $times$ white defendant',
                             'Low education $times$ race frame',
                             'Black defendant $times$ race frame',
                             'White defendant $times$ race frame',
                             'Low education $times$ black defendant $times$ race frame',
                             'Low education $times$ white defendant $times$ race frame',
                             'Strongly Oppose $|$ Somewhat Oppose', 
                             'Somewhat Oppose $|$ Somewhat Favor', 
                             'Somewhat Favor $|$ Strongly Favor'))


##############################################
### Heterogeneous Treatment: Age 18-29     ###
###                 Table S9               ###
##############################################

# binary indicator for age 18-29
WhiteData$Young <- ifelse(WhiteData$Age == 1, 1, 0)

# ordered probit conditioning on young age
YoungOP <- polr(as.factor(SupportDP) ~ Young*Treat1*Treat2,
                data = WhiteData, 
                method = 'probit')

# result
coeftest(YoungOP)

# Marginal effects
MENoPhotoRaceFrame <- glht(YoungOP,
                           linfct = c("Treat21 + Young:Treat21 = 0"))
summary(MENoPhotoRaceFrame)
# power = 0.947
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Young == 1 & WhiteData$TreatmentFactor == 'NP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Young == 1 & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.07142/0.14630,
             sig.level = 0.05,
             alternative = 'two.sided')

MEBlackPhotoRaceFrame <- glht(YoungOP,
                              linfct = c("Treat21 + Young:Treat21 + Treat11:Treat21 + Young:Treat11:Treat21 = 0"))
summary(MEBlackPhotoRaceFrame)
# power = 0.820
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Young == 1 & WhiteData$TreatmentFactor == 'BP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Young == 1 & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.1252/0.1527,
             sig.level = 0.05,
             alternative = 'two.sided')

MEWhitePhotoRaceFrame <- glht(YoungOP,
                              linfct = c("Treat21 + Young:Treat21 + Treat12:Treat21 + Young:Treat12:Treat21 = 0"))
summary(MEWhitePhotoRaceFrame)
# power = 0.988
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Young == 1 & WhiteData$TreatmentFactor == 'WP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Young == 1 & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.08541/0.14822,
             sig.level = 0.05,
             alternative = 'two.sided')

# table output of ordered probit
texreg(YoungOP, stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Age 18-29',
                             'Black defendant',
                             'White defendant',
                             'Race frame',
                             'Age 18-29 $\\times$ black defendant',
                             'Age 18-29 $\\times$ white defendant',
                             'Age 18-29 $\\times$ race frame',
                             'Black defendant $\\times$ race frame',
                             'White defendant $\\times$ race frame',
                             'Age 18-29 $\\times$ black defendant $\\times$ race frame',
                             'Age 18-29 $\\times$ white defendant $\\times$ race frame',
                             'Strongly Oppose $|$ Somewhat Oppose', 
                             'Somewhat Oppose $|$ Somewhat Favor', 
                             'Somewhat Favor $|$ Strongly Favor'))

##############################################
### Heterogeneous Treatment: Age 30-59     ###
###                 Table S10              ###
##############################################

# binary indicator for age 30-59
WhiteData$Middle <- ifelse(WhiteData$Age > 1 & WhiteData$Age < 5, 1, 0)

# ordered robit conditioning on middle age
MiddleOP <- polr(as.factor(SupportDP) ~ Middle*Treat1*Treat2,
                 data = WhiteData, 
                 method = 'probit')

# result
coeftest(MiddleOP)

# Marginal effects
MENoPhotoRaceFrame <- glht(MiddleOP,
                           linfct = c("Treat21 + Middle:Treat21 = 0"))
summary(MENoPhotoRaceFrame)
# power = 0.999
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Middle == 1 & WhiteData$TreatmentFactor == 'NP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Middle == 1 & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.08574/0.12039,
             sig.level = 0.05,
             alternative = 'two.sided')

MEBlackPhotoRaceFrame <- glht(MiddleOP,
                              linfct = c("Treat21 + Middle:Treat21 + Treat11:Treat21 + Middle:Treat11:Treat21 = 0"))
summary(MEBlackPhotoRaceFrame)
# power = 1
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Middle == 1 & WhiteData$TreatmentFactor == 'BP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Middle == 1 & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.1183/0.1223,
             sig.level = 0.05,
             alternative = 'two.sided')

MEWhitePhotoRaceFrame <- glht(MiddleOP,
                              linfct = c("Treat21 + Middle:Treat21 + Treat12:Treat21 + Middle:Treat12:Treat21 = 0"))
summary(MEWhitePhotoRaceFrame)
# power = 0.621
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Middle == 1 & WhiteData$TreatmentFactor == 'WP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Middle == 1 & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.02937/0.11894,
             sig.level = 0.05,
             alternative = 'two.sided')

# table output of ordered probit
texreg(MiddleOP, stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Age 30-59',
                             'Black defendant',
                             'White defendant',
                             'Race frame',
                             'Age 30-59 $\\times$ black defendant',
                             'Age 30-59 $\\times$ white defendant',
                             'Age 30-59 $\\times$ race frame',
                             'Black defendant $\\times$ race frame',
                             'White defendant $\\times$ race frame',
                             'Age 30-59 $\\times$ black defendant $\\times$ race frame',
                             'Age 30-59 $\\times$ white defendant $\\times$ race frame',
                             'Strongly Oppose $|$ Somewhat Oppose', 
                             'Somewhat Oppose $|$ Somewhat Favor', 
                             'Somewhat Favor $|$ Strongly Favor'))


##############################################
### Heterogeneous Treatment: Age 60+       ###
###                 Table S11              ###
##############################################

# binary indicator for age 60+
WhiteData$Old <- ifelse(WhiteData$Age > 4, 1, 0)

# ordered probit conditioning on old age
OldOP <- polr(as.factor(SupportDP) ~ Old*Treat1*Treat2,
              data = WhiteData, 
              method = 'probit')
# result
coeftest(OldOP)

# Marginal effects
MENoPhotoRaceFrame <- glht(OldOP,
                           linfct = c("Treat21 + Old:Treat21 = 0"))
summary(MENoPhotoRaceFrame)
# power = 0.315
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Old == 1 & WhiteData$TreatmentFactor == 'NP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Old == 1 & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.2067/0.4016,
             sig.level = 0.05,
             alternative = 'two.sided')

MEBlackPhotoRaceFrame <- glht(OldOP,
                              linfct = c("Treat21 + Old:Treat21 + Treat11:Treat21 + Old:Treat11:Treat21 = 0"))
summary(MEBlackPhotoRaceFrame)
# power = 0.693
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Old == 1 & WhiteData$TreatmentFactor == 'BP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Old == 1 & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.3741/0.4290,
             sig.level = 0.05,
             alternative = 'two.sided')

MEWhitePhotoRaceFrame <- glht(OldOP,
                              linfct = c("Treat21 + Old:Treat21 + Treat12:Treat21 + Old:Treat12:Treat21 = 0"))
summary(MEWhitePhotoRaceFrame)
# power = 0.951
pwr.t2n.test(n1 = nrow(WhiteData[WhiteData$Old == 1 & WhiteData$TreatmentFactor == 'WP_BF', ]),
             n2 = nrow(WhiteData[WhiteData$Old == 1 & WhiteData$TreatmentFactor == 'NP_NF', ]),
             d = 0.4513/0.3587,
             sig.level = 0.05,
             alternative = 'two.sided')

# table output of ordered probit
texreg(OldOP, stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Age 60+',
                             'Black defendant',
                             'White defendant',
                             'Race frame',
                             'Age 60+ $\\times$ black defendant',
                             'Age 60+ $\\times$ white defendant',
                             'Age 60+ $\\times$ race frame',
                             'Black defendant $\\times$ race frame',
                             'White defendant $\\times$ race frame',
                             'Age 60+ $\\times$ black defendant $\\times$ race frame',
                             'Age 60+ $\\times$ white defendant $\\times$ race frame',
                             'Strongly Oppose $|$ Somewhat Oppose', 
                             'Somewhat Oppose $|$ Somewhat Favor', 
                             'Somewhat Favor $|$ Strongly Favor'))

######################
### Balance Checks ###
###   Table S13    ###
######################

PIDBal <- lm(PID7 ~ TreatmentFactor, data = WhiteData)
IdeolBal <- lm(Ideol ~ TreatmentFactor, data = WhiteData)
EducBal <- lm(Educ ~ TreatmentFactor, data = WhiteData)
GenderBal <- lm(Gender ~ TreatmentFactor, data = WhiteData)
IncomeBal <- lm(Income ~ TreatmentFactor, data = WhiteData)
AgeBal <- lm(Age ~ TreatmentFactor, data = WhiteData)

summary(PIDBal)
summary(IdeolBal)
summary(EducBal)
summary(GenderBal)
summary(IncomeBal)
summary(AgeBal)

